library(Seurat)
library(tidyseurat)
library(tidyverse)
library(harmony)
library(SingleR)
library(tidybulk)
library(ggpubr)
source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bulk.R')

# zuniga 2022 sc---------
zuniga_path <- list.files('mission/exercise/zuniga2022/',full.names = T)

zuniga_mtx <- zuniga_path |>
  str_extract('\\d_.+') |>
  str_replace('_', '-') |>
  set_names(zuniga_path, nm = _) |>
  Read10X(strip.suffix = T)

sobj <- zuniga_mtx |>
  CreateSeuratObject(min.cells = 3, min.features = 200, project = 'Zuniga22exercise')

sobj$mito.ratio <- sobj |>
  PercentageFeatureSet('^MT-')

sobj |>
  VlnPlot('mito.ratio')

sobj <- sobj |>
  quick_process_seurat()

sobj |> count(orig.ident)

monaco <- celldex::MonacoImmuneData()

monaco$label.main |> unique()

sobj <- sobj |>
  mark_cell_type_singler(monaco, fine_label = T, new_label = 'monaco_fine')

sobj |> DimPlot(group.by = 'singler_label')

sobj <- sobj |>
  mutate(group = str_extract(orig.ident, '[:alpha:]+') |>
           case_match('Pre' ~ '1-Pre',
                      'During' ~ '2-During',
                      'Post' ~ '3-Post'))

sobj |>
  mutate(crossed = str_c(singler_label, '_', group)) |>
  DotPlot('FCGR2B', group.by = 'crossed') +
  ylab('Cell type groups')

## B cells
bcells <- sobj |>
  filter(singler_label == 'B cells')

bcells |> VlnPlot('FCGR2B', group.by = 'group') +
  stat_summary(fun = ExpMean, geom = 'crossbar', width = .7) +
  labs(x = 'Group', title = 'FCGR2B expression in B cells')

bcells |>
  SetIdent(value = 'group') |>
  FindAllMarkers(features = 'FCGR2B', logfc.threshold = 0)

bcells |>
  FindMarkers(features = 'FCGR2B',
              ident.1 = 'Post',
              group.by = 'group',
              logfc.threshold = 0)

## DC
sobj |>
  filter(str_detect(singler_label, 'Dend')) |>
  SetIdent(value = 'group') |>
  FindAllMarkers(features = 'FCGR2B', logfc.threshold = 0)

# Hu 2020 parkinson -----------
hu20_meta <- read_delim('mission/exercise/GSE124676_series_matrix.txt.gz', comment = '!Series')

hu20_tidy <- hu20_meta |>
  filter(str_detect(a19_before, 'MET|GSM')) |>
  head(2) |>
  mutate(`!Sample_title` = c('acc','exercise')) |>
  column_to_rownames("!Sample_title") |>
  t() |>
  as.data.frame() |>
  rownames_to_column('time') |>
  as_tibble() |>
  mutate(time = str_extract(time, 'before|after'),
         exercise = str_remove(exercise, 'exercise: '))

## ncbi raw count
urld <- "https://www.ncbi.nlm.nih.gov/geo/download/?format=file&type=rnaseq_counts"
path <- paste(urld, "acc=GSE124676", "file=GSE124676_raw_counts_GRCh38.p13_NCBI.tsv.gz", sep="&");
hu20_count <- read_tsv(path)

## ncbi provided annot
apath <- paste(urld, "type=rnaseq_counts", "file=Human.GRCh38.p13.annot.tsv.gz", sep="&")
annot <- read_tsv(apath)

tdb <- annot[,1:2] |>
  right_join(hu20_count) |>
  pivot_longer(-c(1:2), names_to = 'acc') |>
  left_join(hu20_tidy) |>
  tidybulk(.sample = acc, .transcript = Symbol, .abundance = value)

tdb_tc <- tdb |>
  filter()

tdb_tc <- tdb_tc |>
  preproc_bulk(time)

tdb_tc |>
  plot_qc_bulk(value_scaled, time)

tdb_dea <- tdb_tc |>
   test_differential_abundance(~ 0 + time,
                               contrasts = c('timeafter - timebefore'),
                               omit_contrast_in_colnames = T)

tdb_dea |>
  keep_abundant() |>
  pivot_transcript() |>
  filter(Symbol == 'FCGR2B') |>
  select(Symbol, exercise, logFC, FDR)

tdb_tc |>
  filter(Symbol == 'FCGR2B' & value_scaled < 1000) |>
  mutate(time = fct_relevel(time, 'before')) |>
  ggplot(aes(time, value_scaled, color = time)) +
  geom_boxplot() +
  facet_grid(~exercise) +
  expand_limits(y = 0) +
  theme_pubr() +
  labs(x = 'Time point', y = 'Normalized expression',
       title = 'Parkinson disease patient PBMC FCGR2B expression on exercise')

tdb_tc |>
  filter(Symbol == 'FCGR2B' ) |>
  pull(acc)

tdb_tc |>
  filter(.abundant) |>
  summarise(total = sum(value), .by = acc) |>
  ggplot(aes(total, fct_reorder(acc, total))) +
  geom_point()

## MET ---------
tdb_met <- tdb |>
  filter(exercise == 'MET' & acc != 'GSM3539261')

tdb_met <- tdb_met |>
  preproc_bulk(time)

tdb_met |>
  plot_qc_bulk(value_scaled, time)

tdb_dea <- tdb_met |>
  test_differential_abundance(~ 0 + time,
                              contrasts = c('timeafter - timebefore'),
                              omit_contrast_in_colnames = T)

tdb_dea |>
  keep_abundant() |>
  pivot_transcript() |>
  filter(Symbol == 'FCGR2B') |>
  select(Symbol, exercise, logFC, FDR)

tdb_met |>
  filter(Symbol == 'FCGR2B' & acc != 'GSM3539261') |>
  mutate(time = fct_relevel(time, 'before')) |>
  ggplot(aes(time, value_scaled, color = time)) +
  geom_boxplot() +
  facet_grid(~exercise) +
  theme_pubr() +
  labs(x = 'Time point', y = 'Normalized expression',
       title = 'Parkinson disease patient PBMC FCGR2B expression on exercise')

tdb_met |>
  filter(Symbol == 'FCGR2B') |>
  slice_max(value_scaled, n = 5)

tdb_met |>
  filter(.abundant) |>
  summarise(total = sum(value), .by = acc) |>
  ggplot(aes(total, fct_reorder(acc, total))) +
  geom_point()

# larocca 2022 ----------
laro22 <-
  readxl::read_excel('mission/exercise/GSE206505_VO2-counts-matrix.xlsx')

laro22 <- laro22 |>
  rename(gene = ...1) |>
  pivot_longer(-1) |>
  mutate(
    time = ifelse(
    (str_extract(name,'\\d+') |> as.numeric()) %% 2 != 0,
    'before', 'after'),
    improve = ifelse(str_detect(name, 'NR'), 'no', 'yes'))

tdb.laro <- laro22 |>
  filter(!is.na(value)) |>
  tidybulk(.sample = name,
           .transcript = gene,
           .abundance = value)

tdb.laro <- tdb.laro |>
  preproc_bulk(time)

tdb.laro |>
  plot_qc_bulk(scaled_abundance = value_scaled,
               group = time)

tdb.vo2 <- tdb.laro |>
  test_differential_abundance(~ 0 + time,
                              contrasts = c('timeafter - timebefore'),
                              omit_contrast_in_colnames = T)

tdb.vo2 |>
  pivot_transcript() |>
  filter(gene == 'FCGR2B')

gg.laro <- tdb.laro |>
  filter(gene == 'FCGR2B') |>
  mutate(time = fct_relevel(time, 'before')) |>
  ggplot(aes(time, value_scaled, color = time)) +
  geom_boxplot() +
  theme_pubr() +
  labs(x = 'Time point', y = 'Normalized expression',
       title = 'Healthy women 16 week (n=30)',
       subtitle = 'GSE206505')

## vo2 improved --------
tdb.improv <- laro22 |>
  filter(!is.na(value) & improve != 'yes') |>
  tidybulk(.sample = name,
           .transcript = gene,
           .abundance = value)

tdb.improv <- tdb.improv |>
  preproc_bulk(time) |>
  test_differential_abundance(~ 0 + time,
                              contrasts = c('timeafter - timebefore'),
                              omit_contrast_in_colnames = T)

tdb.improv |>
  pivot_transcript() |>
  filter(gene == 'FCGR2B')

tdb.improv |>
  filter(gene == 'FCGR2B') |>
  mutate(time = fct_relevel(time, 'before')) |>
  ggplot(aes(time, value_scaled, color = time)) +
  geom_boxplot() +
  theme_pubr() +
  labs(x = 'Time point', y = 'Normalized expression',
       title = 'Healthy women PBMC FCGR2B expression on exercise')

## ncbi realigned counts ------
laro22.re <- download_ncbi_counts('GSE206505')

laro22.meta <- read_delim('mission/exercise/larocca_meta.txt',
                          col_names = c('name', 'sample'))

laro22.re <- laro22.re |>
  select(-2) |>
  pivot_longer(-1) |>
  mutate(time = ifelse((str_extract(name, '.$') |> as.numeric()) %% 2 == 0,
                       'after',
                       'before'))

laro22.re <- laro22.re |>
  left_join(laro22.meta) |>
  mutate(improve = str_detect(sample, 'NR', negate = TRUE))

tdb.laro.re <- laro22.re |>
  filter(!improve) |>
  tidybulk(.sample = name,
           .transcript = Symbol,
           .abundance = value)

tdb.laro.re <- tdb.laro.re |>
  preproc_bulk(time)

tdb.laro.re <- tdb.laro.re |>
  test_differential_abundance(~ 0 + time,
                              contrasts = c('timeafter - timebefore'),
                              omit_contrast_in_colnames = T)

tdb.laro.re |>
  pivot_transcript() |>
  filter(Symbol == 'FCGR2B')

tdb.laro.re |>
  filter(Symbol == 'FCGR2B') |>
  mutate(time = fct_relevel(time, 'before')) |>
  ggplot(aes(time, value_scaled, color = time)) +
  geom_boxplot() +
  theme_pubr() +
  labs(x = 'Time point', y = 'Normalized expression',
       title = 'Healthy women PBMC FCGR2B expression on exercise')

# contripois 2020 ----------
cont20 <-
  data.table::fread('mission/exercise/Transcriptomics_VST_excl_3participants.csv')

cont20

cont20.meta <- read_delim('mission/exercise/contripois_meta.txt') |>
  filter(Event == 'Exercise')

cont20.meta <- cont20.meta |>
  select(VisitID, Event_Note1) |>
  dplyr::rename(Sample_ID = VisitID)

cont20.code <- cont20.meta |>
  mutate(suffix = str_extract(Sample_ID, '-.+')) |>
  dplyr::count(suffix, Event_Note1) |>
  select(-n)

cont20.tidy <- cont20 |>
  pivot_longer(-1) |>
  mutate(suffix = str_extract(Sample_ID, '-.+')) |>
  left_join(cont20.code)

cont20.fcgr2b <- cont20.tidy |>
  mutate(ensembl = str_remove(name, '\\..+')) |>
  ensembl_to_symbol(ensembl) |>
  filter(transcript == 'FCGR2B')

cont20.fcgr2b |>
  filter(Event_Note1 %in% c('Baseline', '1 hour')) |>
  mutate(time = fct_relevel(Event_Note1, 'Baseline')) |>
  ggplot(aes(time, value)) +
  geom_boxplot() +
  geom_jitter(width = .1, height = 0)

cont20.ciber <- cont20.tidy |>
  mutate(ensembl = str_remove(name, '\\..+')) |>
  ensembl_to_symbol(ensembl) |>
  filter(Event_Note1 %in% c('Baseline', '1 hour')) |>
  select(Sample_ID, transcript, value) |>
  pivot_wider(id_cols = transcript, names_from = Sample_ID,
              values_from = value, values_fn = sum) |>
  dplyr::rename(Gene = transcript)

cont20.ciber |>
  filter(!is.na(Gene)) |>
  select(1:31) |>
  write_tsv('mission/exercise/cont20.ciber.p1.tsv')

# array data----
## Thompson 2010 ------
library(GEOquery)
library(tidySummarizedExperiment)

gset <- getGEO("GSE12385", AnnotGPL = TRUE) |>
  pluck(1) |>
  makeSummarizedExperimentFromExpressionSet()

id2b <- gset@rowRanges@elementMetadata |>
  as_tibble() |>
  filter(Gene.symbol == 'FCGR2B') |>
  pull(ID)

thomp10.2b <- gset |>
  filter(.feature %in% id2b) |>
  select(title, exprs) 

gg.thmp10 <- thomp10.2b |>
  mutate(time = str_extract(title, '\\d+weeks|basal') |>
           fct_relevel('basal'),
         il6 = str_extract(title, 'low|high')) |>
  filter(time != '26weeks') |>
  ggplot(aes(time, log2(exprs), color = time)) +
  geom_boxplot() +
  theme_pubr() +
  labs(title = 'Healthy men 24 week (n=12)',
       subtitle = 'GSE12385')

## Rampersaud 2013 ---------
gset <- getGEO("GSE34788", AnnotGPL = TRUE) |>
  pluck(1) |>
  makeSummarizedExperimentFromExpressionSet()

id2b <- gset@rowRanges@elementMetadata |>
  as_tibble() |>
  filter(Gene.symbol == 'FCGR2B') |>
  pull(ID)

rampr13.2b <- gset |>
  filter(.feature %in% id2b) |>
  select(title, exprs) 

gg.ramp <- rampr13.2b |>
  mutate(time = str_extract(title, 'BEFORE|AFTER') |>
           fct_relevel('BEFORE')) |>
  ggplot(aes(time, exprs, color = time)) +
  geom_boxplot() +
  theme_pubr() +
  labs(title = 'Healthy women 12 week (n=60)',
       subtitle = 'GSE34788',
       y = 'Normalized expression')

## Dias 2015 ----------
gset <- getGEO("GSE57999", AnnotGPL = TRUE) |>
  pluck(1) |>
  makeSummarizedExperimentFromExpressionSet()

id2b <- gset@rowRanges@elementMetadata |>
  as_tibble() |>
  filter(Gene.symbol == 'FCGR2B') |>
  pull(ID)

dias15.2b <- gset |>
  filter(.feature %in% id2b) |>
  select(title, exprs) 

gg.dias <- dias15.2b |>
  mutate(time = str_extract(title, 'pre|post') |>
           fct_relevel('pre')) |>
  ggplot(aes(time, exprs, color = time)) +
  geom_boxplot() +
  theme_pubr() +
  labs(title = 'Healthy men 18 week (n=13)',
       subtitle = 'GSE57999',
       y = 'Normalized expression')

gg.laro + gg.ramp + gg.thmp10 + gg.dias & NoLegend()
